Role of spatial inhomogeneity in dissociation of trapped molecular condensates 
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We theoretically analyze dissociation of a harmonically trapped Bose-Einstein condensate of 
molecular dimers and examine how the spatial inhomogeneity of the molecular condensate affects 
the conversion dynamics and the atom-atom pair correlations in the short-time limit. Both fermionic 
and bosonic statistics of the constituent atoms are considered. Using the undepleted molecular-field 
approximation, we obtain explicit analytic results for the asymptotic behavior of the second-order 
correlation functions and for the relative number squeezing between the dissociated atoms in one, 
two and three spatial dimensions. Comparison with the numerical results shows that the analytic 
approach employed here captures the main underlying physics and provides useful insights into the 
dynamics of dissociation for conversion efficiencies up to 10%. The results show explicitly how the 
strength of atom-atom correlations and relative number squeezing degrade with the reduction of the 
size of the molecular condensate. 
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I. INTRODUCTION 

Pair production of photons has been the key mecha- 
nism in a series of landmark experiments in quantum op- 
tics. Of particular importance have been the experiments 
with pair correlated photons from an atomic radiative 
cascade [l[ and from parametric down-conversion 0, Q > 
leading to demonstrations of violations of Bell's inequal- 
ities and the Einstein-Podolsky-Rosen (EPR) entangle- 
ment 043 ■ Performed initially with discrete polarization 
states of photons — as a realization of Bohm's version Q 
of the EPR gedankenexperiment @ — the EPR correla- 
tions have been later demonstrated in their original ver- 
sion, that is, for a system of observables with a continuous 
spectrum [10|. In this case, the role of a pair of canon- 
ically conjugate variables is taken by quadrature ampli- 
tudes of spatially separated signal and idler beams gener- 
ated using a nondegenerate parametric down-conversion 
in a cavity. 

In atom optics, one of the simplest mechanisms for 
atom pair production can be realized via dissociation 
of diatomic molecules [n] - [T5| . When performed using 
a molecular Bose-Einstein condensate (BEC) as the ini- 
tial state and assuming that the constituent atoms are 
bosons, the process lends itself [l6l - [l9| as a direct matter- 
wave analog of optical parametric-down conversion of an 
intense laser light in a x -nonlinear crystal. Owing to 
this analogy, one can envisage that future experiments on 
dissociation of molecular BECs can lead to atom optics 
demonstrations of EPR correlations [20l - [22j and related 
tests of Bell's inequalities HHH]. 

The most important ingredients for such experiments 
are: (i) the creation of quantum degenerate samples 
of stable, long-lived molecules, preferably in a rovibra- 
tional ground state (such as in recently demonstrated 
experiments (25l - i27| ). and (ii) the development of tech- 
niques for measuring atom-atom correlations fl3|, l28l - [34| . 
Demonstrations of continuous-variable EPR correlations 



in atomic quadratures will additionally require the mea- 
surement of matter-wave quadrature amplitudes using 
stable, mode-matched local oscillator fields as phase ref- 
erences (2ll [35| . This is a challenging task yet to be 
realized experimentally. In the case of fermionic statis- 
tics of the constituent atoms, dissociation of a molecular 
BEC offers a new paradigm of fermionic quantum atom 
optics [36l [37l | and new opportunities for extensions of 
fundamental tests of quantum mechanics to ensembles of 
neutral fermionic atoms. 

The purpose of this article is a qualitative and quanti- 
tative understanding of the simplest density-density and 
atom number correlations in spontaneous dissociation of 
a molecular BEC in the short-time limit. Strong atom- 
atom correlations and relative number squeezing can be 
regarded as precursors of more complex EPR correla- 
tions. In the short-time limit, the converted fraction 
of molecules into constituent atoms is small and one 
can employ the undepleted molecular-field approxima- 
tion, which was successfully used in theoretical descrip- 
tions of parametric down-conversion in quantum optics 
[HI . Even though a less restrictive description of dissoci- 
ation dynamics can be accomplished via state-of-the-art 
numerical techniques, such as the first-principles simula- 
tions u sing the positive- P and Gaussian representations 
[3^421 , l78fl the truncated Wigner function approach [42| , 
the Hartree-Fock-Bogoliubov method [37L 142 . |43| . and 
a generalized Gross-Pitaevskii (GP) equation [4J|, the 
main motivation of the present work is to obtain ap- 
proximate analytic results which are less computation- 
ally expensive and have the intrinsic appeal of analytic 
simplicity. By comparing these results with the numer- 
ical ones, one can verify their accuracy and the range 
of validity in different parameter regimes. For a related 
problem of atomic four- wave mixing via condensate colli- 
sions |34f. a similar analytic treatment has been employed 
in Ref. [45| , with the results in the short-time limit com- 
paring very well with the first-principles simulations in 
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the positive-P representation (46| . 

The main question that we address here is the influ- 
ence of the spatial inhomogeneity of the source (molec- 
ular) condensate on the strength of atom-atom correla- 
tions. This is an important problem for modeling realis- 
tic experiments performed using trapped, inhomogeneous 
condensates with interactions. For example, for a har- 
monically trapped system, the spatial inhomogeneity and 
the multimode character of the problem enters through 
the shape of the density profile of the initial condensate 
which — in the Thomas-Fermi (TF) limit — is close to an 
inverted parabola. The shape of the source in turn de- 
termines the extent to which the effect of mode-mixing 
affects the spatial structure of atom-atom correlations. 
For strong inhomogeneity, the effect can be quite detri- 
mental and can reduce the strength of correlations sig- 
nificantly. Our analytic results allow us to quantify these 
effects in a relatively simple way, as the short-time limits 
to correlation functions are obtained in a closed explicit 
form in terms of Bessel functions. In a recent rapid com- 
munication [I?}, we presented some of these results for 
dissociation in a one-dimensional (ID) geometry; in the 
present work, we give the details of derivations and ex- 
tend the results to 2D and 3D systems. 

In addition, our comparison between the results for 
bosonic and fermionic atoms allows for a demonstration 
of striking differences in the dissociation dynamics that 
depend inherently on the difference in quantum statis- 
tics. For short durations of dissociation, which produce 
low-density atomic clouds, these differences are not ob- 
structed by the s-wave scattering interactions and have 
been recently studied in the context of directionality ef- 
fects due to Bose enhancement and Pauli blocking in the 
dissociation of elongated molecular condensates [48| . 

Related recent studies of molecular dissociation are 
concerned with the role of confinement on the stability of 
the molecular BEC [4£| , the effect of magnetic- field fluc- 
tuations and modulations on the dissociation dynamics 
near a Feshbach resonance |50l|. the dynamics of disso- 
ciation in optical lattices 1511 . dissociation of molecules 
prepared in a vortex state [52J , connection of dissociation 
with the generators of the SU(1,1) and SU(2) Lie alge- 
bras [53j], loss of atom- molecule coherence due to phase 
diffusion (54|, as well as the use of molecular dissocia- 
tion as a probe of two-body interactions [55[ , collisional 
resonances (5(| , and spectroscopic properties of Feshbach 
resonance molecules |57l . |58| . 

The article is organized as follows. In Sec. |TT] we in- 
troduce the model for dissociation of molecular dimers 
consisting of either two distinguishable bosonic atoms or 
two fermionic atoms in different spin states. In Sec. IIIII 
we formulate the same problem for the case of two indis- 
tinguishable bosonic atoms. In Sec. IIVI we present the 
results of a numerical analysis of molecular dissociation 
in ID within the undepleted molecular-field approxima- 
tion. In Sec. [V] we develop an analytic approach for 
the short-time asymptotic behavior and obtain explicit 
results for the atom-atom pair correlation functions and 



the relative number squeezing in ID, 2D, and 3D ge- 
ometries, for TF and Gaussian density profiles of the 
molecular BEC. Throughout these sections, we discuss 
the validity of the approximate model (with no molec- 
ular depletion) by comparing the relevant results with 
those obtained using first-principles positive-P simula- 
tions for bosons, in which the molecular-field dynamics 
and its depletion is treated quantum mechanically. The 
details of the positive-P method are given in Sec. I VII 
Finally, in the same Sec. I VII we incorporate the effects 
of s-wave scattering interactions and analyze the system 
using the truncated Wigner approach. We conclude the 
article with the summary Sec. IVHI 



II. DISSOCIATION INTO DISTINGUISHABLE 
BOSONIC OR FERMIONIC ATOM PAIRS 

To model the dissociations of a Bose-Einstein con- 
densate of diatomic molecules into pairs of constituent 
atoms, we start with the following effective quantum field 
theory Hamiltonian, in a rotating frame [59j: 

H= [d D x\ ^ T |V$,| 2 + ^A($ t 1 $ 1 +$t$ 2 ) 

J I i=0,l,2 1 

-tfix(*J*i*a-*l*I*o)}. (1) 

Here we assume that the molecules [described by the 
field operator ^o(x,t)] are made of either two distin- 
guishable bosonic atoms or two fermionic atoms in dif- 
ferent spin states. In both cases, ^(x, i) is a bosonic 
field operator satisfying the standard commutation re- 
lation [§o(x,i),§J(x',i)] = <5 D (x-x'), with D = 1,2 
or 3 corresponding to the dimensionality of the system. 
The atomic field operators, \I/i(x, t) (i = 1,2), satisfy 
either bosonic commutation or fermionic anticommuta- 
tion relations, [$j(x, t), $](x', t)} = 5 i j8 D (x. — x') or 

{$i(x, t),$t( x ',f)} = (5 l:( -(5 D (x- x'), depending on the 
underlying statistics. 

The first term in the Hamiltonian ([!) describes the ki- 
netic energy where the atomic masses are mi and ?B2, 
whereas the molecular mass is too = mi + TO2 . For sim- 
plicity, we consider the case of equal atomic masses (same 
isotope atoms), with mi = TO2 = m and mo = 2m. 

The coupling constant \ = Xd is responsible for coher- 
ent conversion of molecules into atom pairs, for example, 
via optical Raman transitions, an rf transition, or a Fes- 
hbach resonance sweep (see, for example, Refs. (59l - [65| 
and (66l |67| for recent reviews); the microscopic expres- 
sions for x m ID, 2D, and 3D can be found in Ref. (43| . 

The detuning A is defined to give the overall energy 
mismatch 2hA between the free two-atom state in the 
dissociation threshold and the bound molecular state (in- 
cluding the relative frequencies of the Raman lasers or 
the frequency of the rf field; for further details, see Refs. 
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[171 143|). Unstable molecules, spontaneously dissociat- 
ing into pairs of constituent atoms, correspond to A < 0, 
with 2fi,|A| being the total dissociation energy. 

The trapping potential for preparing the initial molec- 
ular BEC — with any residual atoms being removed — 
is omitted from the Hamiltonian since we assume that 
once the dissociation is invoked, the trapping potential 
is switched off, so that the dynamics of dissociation is 
taking place in free space (in ID and 2D geometries, we 
assume that the confinement in the eliminated dimen- 
sions is kept on so that the free-space dynamics refers 
only to the relevant dimension under consideration) . We 
assume that the switching on of the atom-molecule cou- 
pling and switching off of the trapping potential is done 
in the regime of a sudden jump [58J. Accordingly, the 
preparation stage is reduced to assuming a certain initial 
state of the molecular BEC in a trap, after which the 
dynamics is governed by the Hamiltonian ([T]). 

In what follows we initially treat the dynamics of dis- 
sociation in the undepleted molecular condensate ap- 
proximation in which the molecules are represented as 
a fixed classical field. The undepleted molecular approx- 
imation is valid for short-enough dissociation times dur- 
ing which the converted fraction of molecules does not 
exceed ~ 10% [H, [H, [H[. In this regime the disso- 
ciation typically produces low-density atomic clouds for 
which the atom-atom s-wave scattering interactions are 
negligible (3?|; hence, their absence from our Hamilto- 
nian. Additionally, the atom-molecule interactions will 
initially appear as an effective spatially dependent de- 
tuning due to the mean-field interaction energy; this can 
be neglected by operating at relatively large absolute val- 
ues of the dissociation detuning |A| so that the total dis- 
sociation energy 2fi,|A| dominates the mean-field energy 
shifts [111 H|| . As means of verifying the regime of va- 
lidity of these approximations, in Sec. I VII we incorporate 
the effects of molecular depletion and s-wave scattering 
interactions using the positive- P representation and the 
truncated Wigner method. As these numerical methods 
are only applicable to bosons, our comparison and conclu- 
sions are restricted to the case of dissociation into bosonic 
atoms. For the case of fermionic atoms, the development 
of stochastic methods that may facilitate similar compar- 
ison in the future are under development [78| . 



A. Heisenberg equations in the undepleted 
molecular condensate approximation 

The undepleted molecular-field approximation is in- 
voked as follows. Assuming that the molecules are 
in a coherent state initially, with the density profile 
po(x) given by the ground-state solution of the stan- 
dard GP equation in a harmonic trap, we replace the 
molecular-field operator by its coherent mean-field am- 
plitude, $ (x,0) -> ($ (x,0)) = * (x,0) = y/ftffi, 
which we assume is real without loss of generality. We 
can next introduce an effective, spatially dependent cou- 



pling, 



,g(x) = xvVo(x), 



(2) 



and write the Heisenberg equations for the atomic fields 
as follows: 



at 
at 



A V 2 - A 

2m 



2m 



$ 1 (x,t)± 5 (x)^(x,t), 



(3) 



*£(x,t)+0(x)*i(x,t). 



Here and hereafter the + and — signs (in general, upper 
and lower signs) are for bosonic and fermionic atoms, 
respectively. 

Transforming to Fourier space, ^ j (x, t) = 
j d D kaj (k, t) exp(ik • x)/(2tt) d / 2 , where the amplitude 
operators a,j (k, t) satisfy commutation or anticom- 
mutation relations, [<2j(k), at(k')] = 5ij5 D (k — k') 

or {oi(k), aj(k')} = ^-^(k — k'), according to the 
underlying statistics, we can rewrite Eqs. ^ as a set of 
linear operator equations: 



d2i(k, t) 
Jt 



iA k ax(k,t)± J -A^ ff(q+k)o|(q,t) 



(4) 



Here 



(2w) D / 2 



-ik-x 



9(x), 



(5) 



is the Fourier transform of the effective coupling g(x), 
and we have defined A^ = hk 2 / (2m) + A, where k = |k|. 

The general structure of solutions following from Eqs. 
(|4"|), with vacuum initial conditions for the atomic fields, 
can be easily understood if we rewrite these operator 
equations in terms of ordinary differential equations for 
all possible second-order moments of the atomic field op- 
erators. By doing so, one can show that the equations 
for the normal and anomalous densities, (aj (k, t)a,j (k', t)) 
and (2i(k, £)o2(k',t)), together with their complex con- 
jugates, form a closed set and develop nonzero popu- 
lations from the delta-function "seed" terms that orig- 
inate from the following identity: (a 3 '(k, £)at(k',£)) = 

5 D (k - k') ± (aj(k,i)a,-(k',*)). The other second-order 

moments, (aj(k, £)a2(k', t)} and (a,j(k, t)aj(k', t)}, also 
form a closed set; however, they never develop nonzero 
populations if the populations were absent initially. Ac- 
cording to this structure, the only nonzero second-order 
moments are the normal and anomalous densities 

» i (k,k , > t) = (at(k,*)a i (k , ,t)> > j = 1,2, (6) 
m 12 (k,k',i) ee (ai(k,t)a 2 (k',£)), (7) 
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with Hi(k, k',i) = ri2(k, k',t), whereas 
(aJ(M)a 2 (k',t)) = and (a,-(k, t)a,-(k', *)) = 0. 
Since the effective Hamiltonian corresponding to Eqs. 
(|3"|) is quadratic in the field operators, any higher-order 
moments or expectation values of products of creation 
and annihilation operators will factorize according 
to Wick's theorem into products of the normal and 
anomalous densities n_j(k, k',t) and mi2(k,k',£). 



III. DISSOCIATION INTO 
INDISTINGUISHABLE BOSONIC ATOMS 

For completeness, we also analyze dissociation of a 
BEC of molecular dimers made of pairs of indistinguish- 
able bosonic atoms in the same spin state. This is de- 
scribed by the following effective Hamiltonian [6(|, in a 
rotating frame: 



H 



/ I _ 777 j 

J I i=0,l 



-i^($t$2_$t2$ o 



(8) 



Here 4'i(x,i) is the atomic field operator, x is the respec- 
tive atom- molecule coupling (43l [63j . mo and mi = m 
are, respectively, the molecular and atomic masses (mo = 
2m), and A is the detuning corresponding to the total 
dissociation energy of 27i|A|. 

The treatment of this system is essentially the same as 
in the previous case of distinguishable atoms, except that 
the field operators ^(x, i) and a 2 (k, t) in Eqs. © and 
(U]) are replaced, respectively, with ^(x, t) and a}(k, t). 
The corresponding Heisenberg equations now read as 



c"Fi(x,t) 
dt 

and 

<iai(k, t) 
dt ~~ 



HV 2 



2m 



$ 1 (x,i)+.g(x)$ t 1 ( Xi t) ; ( 9 ) 



da 

-*A fc Si(k,t) + J . ^ /a g(q+k)o[(q,t), 

(10) 

whereas the nonzero normal and anomalous densities are 
given by 

ni(k,k',i) = (a t 1 (k,t)ai(k',f)), (11) 
ron(k,k',f) = (ai(k,t)ai(k',t)> > (12) 

where we have omitted the atomic spin index for nota- 
tional simplicity. 



IV. NUMERICAL RESULTS AND DISCUSSION 

In our numerical treatment of the problem (present 
section) we only consider a ID system. This is to make 



the problem computationally tractable. However, based 
on the physical understanding that we develop, we ex- 
pect our results to be at least qualitatively valid for 3D 
systems, as is the case for a 2D problem treated recently 
in Ref. [48[. The analytical results of Sec. [V] on the 
other hand, are obtained for ID, 2D and 3D systems. 

From the structure of Eqs. (f?]) we can easily recog- 
nize the role of mode-mixing in the spatially inhomo- 
geneous treatment compared to the case of a uniform 
molecular condensate. In the uniform case, the Fourier 
transform of the effective coupling go is a delta func- 
tion g(k) = (27r)- D / 2 go<5(k), so that the operator ai(k) 
couples to the conjugate of the partner spin component 
at exactly the opposite momentum, a 2 (— k). Therefore, 
the entire set of coupled equations breaks down into 
pairs of equations that couple only the opposite momen- 
tum components of the two atomic fields with different 
spins. Accordingly, only the diagonal and antidiagonal 
terms of the normal and anomalous densities, rij (k, k) 
and mi2(k, — k), develop nonzero populations as the dis- 
sociation proceeds. 

In the present inhomogeneous case, on the other hand, 
the finite width of the effective coupling <jf(k) implies that 
ai(k) couples not only to a 2 (— k) in Eq. (TJJ, but also to 
a range of momenta in the neighborhood of — k, within 
— k±<5k. This is the origin of mode-mixing. The spread in 
<5k determines the width of the pair correlation between 
the atoms in the two opposite spin states that have equal 
but opposite momenta. The width is ultimately related 
to the momentum width of the source molecular conden- 
sate, as we show later in this article. 

At a qualitative level, the finite width of the pair cor- 
relation at opposite momenta can be understood from 
a simple momentum conservation argument. For a 
molecule at rest, the dissociation produces one atom 
in each spin state satisfying ki + k 2 = 0, and there- 
fore k 2 = — ki. From energy conservation 2ft|A| = 
(?i 2 |ki| 2 + fi, 2 |k 2 | 2 )/2m, the absolute momentum of each 
atom is given by fco = |ki| = |k 2 | = y // 2m\A\/h. The 
same momentum conservation holds in the center-of-mass 
frame of the molecule if it has a finite momentum offset 
<5k due to the initial momentum spread of the condensate. 
In the laboratory frame, this center-of-mass momentum 
offset leads to an offset from ±ko in the momenta of dis- 
sociated atoms, ki = ko + <5k/2 and k 2 = — ko + <5k/2. 
This implies that a pair of atom detectors set to detect 
atoms with these momenta will produce a positive pair 
correlation signal and will therefore contribute to the fi- 
nite width of the order of (5k in the density-density cor- 
relation function (68j . 

We now turn to the quantitative analysis of atom-atom 
correlations. Before we proceed, however, we discuss the 
role of mode-mixing on simpler observables — the atomic 
momentum distribution and mode population dynamics 
in the two spin states. 
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Figure 1: (Color online) Molecular BEC density profiles po{x) 
as ground-state solutions of the Gross-Pitaevskii equation in 
a harmonic trap with longitudinal frequencies u/2n = 1, 3/2, 
and 3 Hz represented, respectively, by the curves (1), (2), and 
(3). In the Thomas- Fermi (TF) regime, the corresponding 
TF radii are (1) R TF ^ 250 /im, (2) R TF ^ 167 ^m, and (3) 
Rtf — 83 ^m. The dashed box illustrates a uniform system 
which is size matched with the inhomogeneous system (1). 
Other physical parameters are given in Ref. f6Sj | - 



A. Momentum distribution and mode population 
dynamics 

In a finite quantization volume, the wave vector k 
is discrete and the plane-wave mode annihilation and 
creation operators a^k = o>j 

(MXAfc^AfcyA/y 1 / 2 and 
aj k = a](k,t)(Afc 2; Afc.yAfc 2 ) 1 / 2 (where j = 1,2 and 
A/c Xi y iZ are the lattice spacings in x, y, and z directions) 
may be organized into a vector a. The Heisenberg equa- 
tions (|4]) may then be written in vector-matrix form as 
da/dt = Ma, where M is a square matrix of complex 
numbers of dimension equal to twice the total number 
of lattice points. The solutions of these linear operator 
equations can be found by numerically computing the 
matrix exponential exp(Mt). The task is relatively sim- 
ple in ID, which is the case that we present here. 

In our numerical analysis, we consider three typical 
examples of the density profiles po(x) of the molecular 
BEC, corresponding to relatively weak, intermediate, and 
strong inhomogeneity. These are shown in Fig. Q] and 
correspond to having different frequencies of the longi- 
tudinal trapping potential along x and the same peak 
density po = po(0). The density profiles shown are given 
by the ground-state solution of the ID GP equation in a 
harmonic trap and can be closely approximated by a TF 
inverted parabola po(x) = po(l — x 2 /R^ F ) for |x| < Rtf 
and po(x) = elsewhere; the three examples shown cor- 
respond to the TF radii of Rtf = 250 pm, Rtf = 167 
pm, and Rtf = 83 pro.. 

To allow for a comparison of the present nonuniform 
treatment with the known analytic solutions of a uniform 
model j36[, we also show a uniform box system (dashed 
line in Fig. [1]), which is size-matched with the largest 
trapped system, curve (1). We choose the size-matched 
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Figure 2: (Color online) (a) Total fractional number of atoms 
Nj(t)/N (0) [Ni(t) = N 2 (t)] as a function of time t/t in 
one of the spin components in dissociation into fermionic or 
bosonic atoms. The different curves (1), (2), and (3) cor- 
respond to dissociation of molecular BECs having density 
profiles as indicated in Fig. [1] with the dashed line refer- 
ring to the uniform case. The results are obtained within 
the undepleted molecular approximation. The dimension- 
less detuning S = Ato is 5 = —9, where the time scale is 
to = 1/xVPo- For X = 3.15 x 10 -2 m 1/2 /s and p = 4 X 10 7 
m _1 , the time scale would be to = 5 ms. Accordingly, the 
absolute detuning is |A| = 1800 s _1 and the resonant mo- 
mentum is ko = y/2m\A\/h ~ 1.5 x 10 6 m -1 in the present 
example [69l |. (b) Total fractional numbers of bosonic atoms 
Nj(t)/No(0) from numerical simulations that take into ac- 
count molecular depletion. The different curves are as in (a), 
with the curve (1) from the undepleted molecular approxi- 
mation shown for comparison. The solid curves are from the 
truncated Wigner method, while the circles are from the ex- 
act positive-P method (see Sec. IVI[) . which in this example 
is limited to simulation duration of t/to ~ 6.5 due to the 
growing sampling errors. 

uniform box to have the same uniform density p u as the 
peak density po of the nonuniform system and the same 
total initial number of molecules A^o- For a simple TF 
parabola, Nq is given by Nq = 4poRtf/3, while for the 
box system Nq = p u L, and therefore we require that the 
box length is L = 4Rtf/3- 

In Fig. HJa) we plot the total number of atoms in each 
spin component, 

Nj(t) = / dknj{k,t), (13) 



(i 




where 



n,j(k,t) =n,j(k,k,t) = (a' ! j (k,t)a j (k,t)) (14) 

is the density distribution in momentum space. The three 
curves (1), (2), and (3) referring to fermionic atoms cor- 
respond, respectively, to dissociation of molecular con- 
densates with the density profiles shown Fig. Q] Plot- 
ted are the fractional atom numbers, Nj(t)/N o (0) [with 
iVi(i) = N2(t)], where A^o(0) in each case is the respec- 
tive total initial number of molecules. For comparison, 
we also show the result for bosonic atoms correspond- 
ing to the molecular density profile (1). The two dashed 
curves are the respective (fermionic or bosonic) results 
for a uniform system, size matched with the molecular 
condensate profile (1). 

The differences between the curves (1), (2), and (3) 
demonstrate the strong dependence of the dissociation 
dynamics on the inhomogeneity of the initial molecular 
condensate. We note that all fermionic examples are in 
the parameter regime where the dynamics is dominated 
by Pauli blocking of individual atomic modes rather than 
by molecular depletion jH}. Accordingly, only a small 
fraction of molecules is converted into atoms and this 
justifies the use of the undepleted molecular approxima- 
tion for longer time scales than in the respective bosonic 
systems. 

In Fig. [U(b) we show the results of numerical simula- 
tions that go beyond the undepleted molecular approxi- 
mation. The simulations are performed for dissociation 
into bosonic atoms using the positive-P representation 
and the truncated Wigner methods (see Sec. IVII for fur- 
ther details). They take into account the conversion dy- 
namics and quantum fluctuations of the molecular field, 
in contrast to the results of Fig. [U[a). As we see, the 
results of the undepleted molecular approximation are 
in excellent agreement with the exact positive-P results 
for dissociation durations up to t/to ~ 3.5 corresponding 
to more than 10% conversion. At later times, the ex- 
act results show the slowing down of the atom number 
growth due to the depletion of the molecular condensate, 
followed by the reverse process of atom-atom recombina- 
tion into molecules. 

To further illustrate the differences between the uni- 
form and nonuniform results, we plot in Fig. [3] the mo- 
mentum distribution of the dissociated fermionic atoms, 
n,j(k,t), as a function of time. The distribution is sym- 
metric around the origin and has two peaks centered 
around the resonant momenta k = ±fco; in Fig. Owe 
only show the spectrum around k = ko- Qualitatively, 
the momentum distribution is similar to the one obtained 
within the uniform treatment except that the oscillation 
of the resonant momentum is no longer periodic and the 
minima do not reach zero [111 |43[ . We recall that in the 
uniform case the oscillations of different plane- wave mode 
occupancies are periodic and are given by (36| 
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Figure 3: (Color online) Atomic momentum distribution 
n,j(k,t) \n\(k,t) = n2{k,t)\ in one of the spin components 
for dissociation into fermionic atoms, as a function of a scaled 
time t/to. The example is for the molecular BEC profile cor- 
responding to the case (1) in Fig. [1] The dimensionless de- 
tuning is S = —9, as in Fig. [2] 
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Figure 4: (Color online) (a) Peak density rij(ko, i) [ni(fco, t) = 
ri2(ko,t)] as a function of time t/to in one of the spin com- 
ponents in dissociation into fermionic or bosonic atoms. The 
three different curves for fermions correspond to the molecular 
BEC profiles (1), (2), and (3) of Fig. \T\ while the oscillating 
dashed curve is the peak atomic density for a uniform system 
that is size matched with (1). The bosonic curve together 
with the respective dashed curve for a size-matched uniform 
system is for the molecular BEC profile (1) of Fig. [1] The 
dimensionless detuning is the same as in Fig. [2] S = —9. (b) 
Comparison of the average mode occupancy in the uniform 
fermionic system (dashed curve) and the respective nonuni- 
form systems (solid curves), corresponding to cases (1), (2), 
and (3) of Fig. [1] The latter three curves are almost indis- 
tinguishable from each other on this scale. 



where go = g(0) is the effective coupling and = 
Hk 2 /(2m) + A = h(k 2 — k 2 )/2m. In the bosonic case, 
we find that the density distributions are closer to the 
respective uniform results of Ref. jH, |3^|, rij.fe(t) = 

and we do not show 



Al t), 



n hk {t) 



9 2 o 



A? 



sin'U/^ + A^ t ), 



(15) 



[g 2 /(gi-Al)]smh 2 (^gJ 
them here. 

In Fig. QJa) we monitor the atomic density at the 
resonant momentum fco as a function of time; the three 
fermionic curves correspond, respectively, to the molecu- 
lar density profiles (1), (2), and (3) of Fig. [TJ whereas the 
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oscillatory dashed line corresponds to the analytic solu- 
tion for a size-matched uniform system (1). The analytic 
result represents the density distribution n,j(ko,t) and is 
obtained from the average occupation number of the reso- 
nant mode, %,fe (t) = sin 2 (got), by converting it into the 
density distribution via i%j(ko,t) = rij t k (t)/ Ak, where 
Afc = 2?t/L is the mode spacing of the size-matched uni- 
form system of length L. The respective bosonic results 
for the largest molecular BEC are also shown for com- 
parison; in this case the uniform analytic result for the 
resonant mode is given by rij_k (t) — sinh 2 (got) (3(1 139| 
and grows exponentially with time due to Bose stimula- 
tion. 

An alternative way of comparing the results for the 
uniform and nonuniform systems is to define the physical 
"modes" of the nonuniform system and to compare their 
average occupation numbers with those obtained in the 
uniform finite size box. In the later case, the natural 
modes of the system are the plane-wave modes which 
coincide with our computational lattice modes. In the 
case of dissociation into fermionic atoms, the plane-wave 
mode occupation numbers are bound to be no more than 
1 by the Pauli exclusion principle [3(| . In the nonuniform 
system, on the other hand, the plane- wave modes are not 
the natural modes of the system and therefore care should 
be taken when defining the physical modes and discussing 
the Pauli exclusion principle. Following Glauber's th eory 
of optical coherence in the context of matter waves [7(| , 
we define the atomic "mode" in the nonuniform system 
using the first-order coherence length, Ak^ coh \ which in 
turn is defined via the first-order correlation function 



Afc( coh ' as the half-width at half maximum becomes less 
appropriate since this definition ignores the correlation 
peaks at large \k — k'\. 

Even though we do not present explicit numerical re- 
sults for the first-order correlation function itself, we 
point out that in the present model the absolute value 
of gfj(k, k', t) is related to the second-order correlation 
function g--(k,k,t) via 



$\k,k\t)\ = yj\g%\ k ,k>,t) 



(17) 



We now turn to the analysis of the second-order corre- 
lation functions for pairs of atoms in the same spin state, 

(2) 

g)j (k,k ,t) (j = 1,2), and in the opposite spin states, 
g[f(k,k',t). 

B. Atom-atom correlations 

1. Distinguishable fermionic or bosonic atoms 

Since the dissociation of diatomic molecules produces 
pairs of atoms in two different spin states which fly apart 
in opposite directions according to the momentum con- 
servation, we expect strong correlation signal for atom 
pairs with nearly equal but opposite momenta. We refer 
to this type of correlation as back-to-back (BB) correla- 
tion and quantify it via Glauber's second-order correla- 
tion function. 



(at(fc,f)ai(fc',*)> 



V^WhWJ) ^ - /a\(k,t)Mk,t))(aUk\t)a 2 (k',t)) 



(18) 



More specifically, we define Afc( coh ' as the distance \k— k'\ 
over which |aj* (fc, k' , t)\ reduces by a factor of 2 from its 
peak value at k — k! = 0. The average mode occupa- 
tion is then given by rij^(t) = nj(k,t)Ak^ coh \ In Fig. 
|D(b) we plot rij t k (t)-, defined in this way, and compare 
it with the occupancies of the plane-wave modes of the 
uniform system (dashed line). As we see, the mode oc- 
cupancy in the nonuniform system deviates substantially 
from the uniform result and remains below the maximum 
occupancy of 1. We note that the nonuniform results 
corresponding to cases (1), (2), and (3) of Fig. Q] al- 
most coincide with each other and the respective three 
curves in Fig. 0Jb) are almost indistinguishable until the 
first oscillation maximum. This implies that the mode 
occupation dynamics depends only on the shape of the 
molecular condensate and not on its size, at least in the 
short-time limit. [For example, for a Gaussian shape of 
the molecular condensate, we find that the mode occu- 
pation dynamics is slightly different.] In the longer time 
limit, the first-order correlation function in the fermionic 
case develops complicated multipeak structure (similar 
to the one seen in the momentum distribution of Fig. 
[3]). In this case, defining the first-order coherence length 



Apart from the normal ordering of the creation and 
annihilation operators, the pair correlation function 

( 2) 

g[ 2 (k,k',t) describes the density-density correlation be- 
tween the momentum components k and k' of pairs of 
atoms in the two spin states. The normalization with re- 
spect to the product of individual densities rij(k,t) = 
(aj(k,t)aj(k,t)) ensures that for uncorrelated states 

( 2) 

g\ 2 (k,k',t) = 1. Due to obvious symmetry considera- 

tions, g[l ) {k,k',t)=g^(k,k',t). 

The second type of correlation expected to be present 
in the system is between pairs of atoms in the same spin 
state, propagating with the nearly same momenta, k ~ 
k! . This type of correlation, which we refer to as collinear 
(CL) correlation, is due to quantum statistical effects and 
represents a manifestation of the Hanbury Brown, and 
Twiss (HBT) effect [H, HI Il3- The CL correlations 
are quantified via the following second-order correlation 
function: 



0$ (*,*',*) 



(a] (k, t)a] (k', t)a,j (k' , t)aj (k, t)) 
(a] (k, t)a 3 (k, t)> (at (k',t)a (k>, t)) 



(19) 
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Figure 5: (Color online) Atom-atom pair correlations at 
t/to = 1 for the case of the molecular BEC profile (f) in Fig. 
[l]and dimensionless detuning 5 — —9. (a) g^ (k, ko, to) as a 
function of k for fermionic and bosonic atom pairs in two dif- 
ferent spin states, showing the peak corresponding to strong 
back-to-back correlation around k — —ko- (b) Pair correla- 
tion for the same-spin atoms (k, ko,to) in dissociation into 
distinguishable (fermionic or bosonic) atom pairs, showing a 
peak (bosons) or a dip (fermions) at k — ko, which corre- 
spond, respectively, to bosonic bunching due to Bose stimula- 
tion or fermionic antibunching due to Pauli blocking, (c) Pair 
correlation g[ i (k, ko, to) in dissociation into indistinguishable 
(same spin state) bosonic atom pairs, showing both back-to- 
back and CL correlation signals at k — — ko and k — ko- 



The linearity of Eqs. ((4]) ensures that one can ap- 
ply Wick's theorem to Eqs. (fP8]l and (fT9|) and factorize 
the fourth-order operator moments into the sum of prod- 
ucts of second-order moments. Noting in addition that 
(a\(k,t)a 2 (k',t)) = and (a j (k,t)a j (k' ,*)) = in the 
present model, we obtain the following results for the 
BB and CL correlations: 

®^ = ^%WJY (21) 

where the + and — signs stand for bosonic and fermionic 
atoms, respectively. 

(2) 

The BB pair correlation g\ 2 (k,ko,t = to) in which 



the momentum of one of the atomic spin components is 
fixed to the resonant momentum ko, while the momen- 
tum k of the opposite spin component is being varied is 
plotted in Fig. (Ha). The two curves correspond to the 
fermionic and bosonic atom statistics, as indicated by the 
labels. In both cases, we see a clear correlation peak at 
k = — ko corresponding to atom pairs with equal but op- 
posite momenta (ko, —ko)- The width of the correlation 
is discussed in Sec. IIV CI 

The major quantitative difference between the present 
nonuniform result and that of a uniform system is that 
the peak value of the BB correlation g\ 2 (— ko, ko, t) be- 
comes smaller than in the uniform system and that the 
correlation function acquires a finite width. As has been 
shown in Refs. (3rl |43|. the strength of the pair correla- 
tion in the uniform system corresponds to its maximum 
possible value for a given occupancy ni yko (t) = n 2tkg (t). 
More specifically, the pair correlation upper bound is 
given by g 12 (— ko, ko, t) = l/ni k (t) in the case of 

(2) 

fermionic atoms and g\ 2 {— ko, ko, t) = 2 + l/ni tko (t) 

(2) 

in the bosonic case, whereas g\ 2 (k, ko, t) = 1 for any 
k =/= —ko in both cases. In other words, the pair cor- 
relation in the uniform system is a Kronecker-like delta 
function, whereas in the nonuniform case it acquires a 
finite width and the peak value is reduced. 

Given that n,j^ (t) in the uniform system is an 
oscillatory function for fermions and can reach zero 
values at certain times, the respective BB corre- 

(2) 

lation g 12 (— ko, ko, t) = l/ni yko (t) is discontinuous 

(2) 

o, k , t) -¥ co] whenever n jyko (t) ->• 0. In contrast 
to this unphysical result, the peak value of the BB corre- 
lation in the nonuniform system is always continuous and 
is shown in Fig. [5] The oscillatorylike structure of the BB 
correlation for fermions comes from the oscillatory behav- 
ior of the atomic density rij(k, t) and the absolute value of 
the anomalous density \mi 2 (ko,~ko,t)\. The oscillations 
in \mi 2 (ko, —kp, t)\ resemble those for a uniform system 
[H|, with \m ka - ko (t)\ 2 = n jyko (t)[l - n jM (t)}, except 
that here they do not reach the minimum and maximum 
values corresponding to perfect harmonic oscillations in 
n jM The BB correlation for the case of bosonic atoms 
is also shown on the same figure (solid red curve) and does 
not display any oscillations. We note that the bosonic 
curve should have been stopped at t/to ~ 3 as the un- 
depleted molecular approximation breaks down beyond 
t/to ^ 3 (the total number of atoms produced beyond 
this point in time corresponds to a conversion of more 
than 10% of the initial number of molecules) , while this is 
not the case for fermionic atoms. We extend the bosonic 
curve to t/to ~ 5 in order to make visible its discrepancy 
with the numerical results based on the exact positive- 
P representation (red circles) and the truncated Wigner 
method (red dash-dotted curve). 

(2) 

^(fcjfco,*) [with 



The CL correlation function 
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Figure 6: (Color online) Back-to-back pair correlations 
(—ko, ko,t) at equal but opposite peak momenta as a 
function of time t/to for dissociation into fermionic (blue 
solid curve) and bosonic (red solid curve) atoms in two spin 
states from the numerical simulations within the undepleted 
molecular-field approximation. The results are for the molec- 
ular density profiles corresponding to case (1) in Fig. [T] The 
dimensionless detuning is S = —9 as in Fig. [5] The dashed 
(black) curve is the short-time analytic result of Eq. (|38[) from 
Sec. [V] the curve is extended to time duration of t/to ~ 1, 
which, strictly speaking, is beyond the expected validity of the 
perturbative theory, t/to "C 1, yet we see a reasonably good 
agreement with the numerical results for t/to ~ 1. The dash- 
dotted (red) curve is from the truncated Wigner approach, 
while the red circles are from the exact positive-P method; 
both methods take into account the molecular-field depletion 
(see Sec. I VII for details). As we see the result from the unde- 
pleted molecular-field approximation (solid red curve, bosons) 
are in excellent agreement with the positive-P and truncated 
Wigner results for durations up to t/to ~ 5. 

atoms we see the expected HBT peak at k = kg, 
with g\ x (fco, ho, t) = 2 due to bosonic stimulation. In 
the fermionic case, on the other hand, we see a dip 

(2) 

fill (ko, fco, t) = corresponding to fermionic antibunch- 
ing due to Pauli blocking. The physics behind these 
effects is the same as the CL correlations between the 
s-wave scattered atoms in a collision of two conden- 
sates of metastable 4 He* atoms as observed in Ref. (34| 
and the local correlations in a thermal cloud of bosonic 
4 He* atoms and an ultracold cloud of fermionic 3 He* 
atoms as observed in Refs. [32| (see also [HI). 

2. Indistinguishable bosonic atoms 

In the case of dissociation of molecules made of bosonic 
atom pairs in the same spin state, the atom-atom pair 
correlation function is given by 

(2)/, ,/ s = (a{ (k, t)a\ (fc', tjfai (k' , tjai (fc, t)) 
' (ol(A,t)oi(A; > t))(ol(fc',*)oi(*' > t)) 
\n 1 (k,k',t)\ 2 + \m 11 (k,k',t)\ 2 

n x {k,t)n r {k' ,t) ' 1 ' 
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Figure 7: (Color online) (a) Back-to-back pair correlation 
functions, g^(k,ko,to), at time t — to for fermions (blue 
solid curve) and bosons (red solid curve), corresponding to 
the molecular BEC profile (f ) of Fig. [T] The magnitudes of 
the correlation functions are scaled to f in order to simplify 
the comparison of the correlation widths and the width of 
the source (dashed curve), that is, the momentum distribu- 
tion of the molecular BEC. (b) CL pair correlation functions, 
9jf(k, fco, to), for fermions (blue solid curve) and bosons (red 
solid curve), for the same parameters as in (a). The fermionic 
curve for g^fk, fco, to) — 1 is inverted to allow for the com- 
parison of the correlation widths. The dashed curve is the 
momentum distribution of the source. 

Here, the BB correlation signal comes from the anoma- 
lous density term |mn(fc, fc', t)\ 2 as it is nonzero only 
for pairs of momenta fc and fc' that are nearly opposite, 
whereas the normal density term |ni(fc, fc', t)\ 2 is vanish- 
ingly small for opposite pairs of momenta. On the other 
hand, the CL correlation signal comes from the normal 
density term |ni(fc, fc', t)\ 2 which is nonzero for pairs of 
nearby momenta, while the anomalous density is vanish- 
ingly small when fc ~ fc'. 

In Fig. (5fc) we plot the pair correlation g^(k, fco, to) 
as a function of fc, and we see the simultaneous presence 
of two peaks: one at fc = —fco representing the BB corre- 
lation due to the momentum conserving pair-production 
process and the second peak at fc = fco corresponding to 
the HBT effect. 



C. Width of the correlation functions 

An important observable in the experiments on atom- 
atom correlations is the width of the correlation func- 
tions. Here we discuss the BB and CL correlation widths 
during relatively short durations of dissociation, cor- 
responding to the range of validity of the undepleted 
molecular-field approximation. 

In Fig. [7| we show the numerical results for the BB 
and CL correlations functions at time t = to and com- 
pare them with the shape of the source molecular BEC in 
momentum space. The results correspond to the molec- 
ular BEC density profile (1) of Fig. [TJ We see that in 
the short-time limit the width and the overall shape of 
the correlation functions for both fermionic and bosonic 
atoms is determined essentially by the width and the 
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Figure 8: (Color online) Widths (solid lines: bosons, red; 
fermions, blue) of the BB and CL pair correlations, tu' BB ' 
and m/ cl ' , relative to the width of the source condensate as a 
function of time, for the molecular BEC profile (1) of Fig. [JJ 



shape of the momentum distribution of the source molec- 
ular BEC. 

For quantitative purposes we define the widths as the 
half-width at half maximum and denote them u/ BB ) and 
W (CL) £ or back-to-back and for CL correlations, re- 
spectively. The width of the source molecular BEC in 
momentum space is denoted as w. The widths of the 
BB and CL correlations relative to the source width as a 
function of time are shown in Fig. [8j The source conden- 
sate here corresponds to the molecular BEC profile (1) 
of Fig. [TJ having the momentum width of w ~ 1.62/i?TF, 
where Rtf is the respective TF radius in x space. We 
see that the BB correlation width starts from the asymp- 
totic value of u/ BB ) = w, corresponding to the limit 
t — > 0, and grows slowly for bosons and decreases for 
fermions. A similar behavior is seen in the CL correla- 
tion width, except that it starts from a higher asymp- 
totic value, w( CL) ~ 1.12w, than the BB width. The 
short-time asymptotic behavior is in agreement with the 
analytic results obtained in Sec. [Vj 

For longer durations of dissociation — while still within 
the range of validity of the undepleted molecular-field 
approximation — we see that the BB correlation width 
grows faster (for bosons) than the CL correlation width, 
indicating the possibility that it may eventually become 
broader than the CL width. A similar behavior has been 
recently found in first-principle simulations of a related 
system of atomic four-wave mixing via condensate col- 
lisions p5j . Moreover, it has been shown that in the 
long-time limit the BB correlation width indeed becomes 
broader than the CL correlation, in agreement with the 
experimental measurements of Ref. [34| . 

For comparison, we have also performed an analysis of 
the dynamics of the BB and CL correlation widths for 
a Gaussian density profile of the source molecular BEC. 
The numerical results are in good agreement with the 
analytic results presented in Sec. IV Bl 



D. Relative number squeezing 

As an alternative measure of the strength of atom- 
atom correlations, we now calculate the variance of rela- 
tive atom number fluctuations for atoms in different spin 
states and with equal but opposite momenta ±fco, 



V ko .-k a {t) = ([A(n 1M -n 2i ^ )] 2 )/A SN , 



(23) 



where Asn is the shot-noise level that originates from 
uncorrelated states. The atom number operators are 
defined by nj,±fc (i) = J K dkrLj(k,t) [with rij(k,i) = 
a,j(k, t)a,j(k, t)], where K is the counting length around 
±fcrj. On a computational lattice the simplest choice of 
K that does not require explicit binning of the signal is 
K = Ak, where Ak is the lattice spacing, and therefore 
rij,±k (t) = rij(±ko,t)Ak. 

The shot-noise level Asn is different for bosons and 
fermions. For the bosonic case, Asn is given by the 
sum of variances of the individual mode occupancies with 
Poissonian statistics (as in the coherent state), implying 
that Asn = (^i,fc ) + (^2,-0- For the fermionic case, 
the sum of the variances of two uncorrelated modes is 

A S n = ("t.fcoX 1 - («i,fe ))+ ("2,-feoX 1 - ("2,-fco)) & 
which is independent of the states of individual modes. 
The relative number variance f2"3"]) in both cases can be 
combined into the following expression: 



V ko ,- ko (t) = 1 - 



x[9i2( k o,-ko,t) -9n (ko,k ,t) 



Afcm(fco, t) 

1 — sAkni(ko, t) 

(2)/ 



(24) 



where s — and s = 1 for bosons and for fermions, re- 
spectively, and we have taken into account that (ni_k ) = 

^ (2) (2) 

("2,-0 and (fco, ko,t) = #22 (~ fc o,-fco,0- Note that 

(2) 

g\i (ko, ko,t) = in the fermionic case due to the Pauli 
exclusion principle. Variance VO-fc W < 1 implies 
squeezing of relative number fluctuations below the shot- 
noise level. 

In Fig. [9] we plot the relative number variance for 
dissociation into fermionic (solid oscillatory curves) and 
bosonic (solid curves that initially follow the dashed hor- 
izontal lines) atom pairs for the three different sizes of 
the molecular BEC density profile of Fig. [TJ The hor- 
izontal dashed lines are the respective analytic results 
(see Sec. [Vj in the short-time limit showing a common 
asymptotic behavior for bosons and fermions in the limit 
t — > 0. In this limit the atomic mode populations are 
much smaller than 1 and the quantum statistical effects 
do not show up. The numerical results for bosonic atoms 
are stopped around t/tg ~ 3.5; at this time the total num- 
ber of atoms produced corresponds to more than 10% 
conversion and the undepleted molecular approximation 
is no longer valid. For the fermionic atoms the converted 
fraction remains less than 10% for the entire time win- 
dow plotted. The circles that follow the bosonic solid 
curves for up to t/to ~ 3.5 are from the exact positive- 
P method (see Sec. IVI[) and demonstrate the validity of 
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Figure 9: (Color online) Relative number variance for dis- 
sociation into fermionic and bosonic atom pairs (full curves) 
within the undepleted molecular approximation. The dashed 
horizontal lines are the analytic results in the short-time 
asymptotic limit [Eq. ([40 [l ] from Sec. [V] The three different 
sets of curves correspond to three different molecular BEC 
profiles of Fig. [T] The circles are from the exact positive- P 
simulations for bosons that take into account the molecular 
depletion. 



the undepleted molecular approximation within this time 
window. 

We see that both bosonic and fermionic cases display 
relative number squeezing that depends strongly on the 
size of the molecular BEC. The larger the molecular BEC 
is, the stronger is the relative number squeezing, imply- 
ing that the degrading role of mode-mixing due to the 
source inhomogeneity is a weaker effect when the molec- 
ular BEC density profile is closer to uniform. The reduc- 
tion of squeezing and the nontrivial oscillatory behavior 
in the fermionic cases reflect the oscillatory behavior of 
the fermionic mode populations and the anomalous den- 
sity \mi2(ka, — fco, t)\, as explained in the discussion of 
the oscillations in the BB peak correlation in Fig. [5J at 
time instances when the populations become close to 
or 1 , the fermionic shot- noise itself becomes vanishingly 
small and therefore the degree of squeezing below the 
shot-noise level diminishes [36| . 

In all three examples of Fig. [9] the lattice spacing is 



chosen as Afc = 1875 



ir/2R$l 



where Rjtp is 



the TF radius of the molecular BEC profile (3) of Fig. 
[TJ This is the largest spacing that is capable of resolving 
the relevant details and correlation lengths in momentum 
space for all three cases, and we see that the degree of 
squeezing is not particularly large in all examples. The 
strongest squeezing is for the case of the largest molecular 
BEC [curve (1)] and is only ~ 14% at short times. 

This relatively modest degrees of squeezing in the ex- 
amples of Fig. [H] can be contrasted to the ideal situation 
of 100% squeezing achievable in a completely uniform 
system that permits analytic solution in the undepleted 
molecular approximation (36[. The strong departure in 
the actual degree of squeezing from the prediction of the 



idealized uniform model is perhaps the strongest mani- 
festation of the role of mode-mixing in realistic inhomo- 
geneous systems. As shown previously [10, [HI, squeezing 
can be enhanced by binning the atomic signal into bins of 
larger size, in which case Eq. (|24p is no longer applicable. 



ANALYTIC TREATMENT IN THE 
SHORT-TIME LIMIT 



In this section we present the results of an analytic 
treatment of the problem of molecular dissociation in 
the short-time limit in ID, 2D, and 3D. Our approach 
is based on perturbative expansion in time, starting 
with the operator equations of motion in the undepleted 
molecular-field approximation [Eq. Q]. Even though 
the present results are applicable for even shorter time 
scales than those of the numerical treatment of the pre- 
vious section, their analytic transparency provides useful 
insights into the problem of atom-atom correlations as it 
has recently been demonstrated for a closely related prob- 
lem of atomic four- wave mixing via condensate collisions 
[45| . For molecular dissociation, the present analytic ap- 
proach in nonuniform ID systems has been employed in 
Ref. here we present the details of derivations and 

extend the results to 2D and 3D systems. 

The short-time perturbative treatment is based on the 
Taylor expansion in time, up to terms of order t 2 , 



aj(k,i) = Oj(k,0)-M 



daj(k, t) 



dt 



t 2 d 2 aj(k,t) 



dt 2 



(25) 

which is valid for t -C to, where to = 1/x^/Po is the 
characteristic time scale. We recall that, in the defi- 
nition of to for ID, 2D, and 3D systems, the coupling 
constant x anc i the molecular BEC peak density po are 
to be understood as their ID, 2D, and 3D counterparts 
(see, e.g., [43|), but we suppress the respective indices 
for simplicity. In all cases, the units of \ an d Po are 
such that l/x\/7>o has units of time. With the pre- 
ceding expansion, one can check that the commutation 
(for bosons) and anticommutation (for fermions) rela- 
tions for the creation and annihilation operators are given 
by [a;(k, t), a](k', t)] T ~ Sij5(k — k'), up to terms of the 
order of t 2 . 

Using the right-hand sides of the generating equations 
of motion ((4]) for calculating the derivative terms in Eq. 
(f2"S")) . we obtain the following expressions for the anoma- 
lous and normal densities in the lowest order in t, in 
D = 1,2, and 3 dimensions: 

|mi 2 (k,k',i)| = |(ai(k,f)a a (k',t)}| 
t 



|p(k + k')| 



(2tt) d / 2 
J ^xe l(k+k ')- x g(x) 



(26) 
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n,(k,k',t) = (at(k,f)aj(k',i)) 



(2 



^yj J ^qg*(q + k)g(q + k') 



(27) 



Recalling Eqs. (J20j) and ([21]) and viewing Eqs. ([26]) 
and ([2"T]) in the context of BB and CL correlation func- 
tions, we see that the width of the CL correlation be- 
tween the same-spin atoms [Eq. (|21[) ] is determined by 
the square of the Fourier transform of the square of the 
effective coupling <?(x). The width of the BB correlation 
[Eq. ([2"0"]) [ between the different spin-state atoms, on the 
other hand, is determined by the square of the Fourier 
transform of g(x). Since the function <?(x) 2 is narrower 
than g(x) and the converse is true for their respective 
Fourier transforms, we immediately deduce that the CL 
momentum correlation is generally broader than the BB 
correlation. These conclusions are true for any shape of 
the source condensate and apply to both bosonic and 
fermionic statistics in the short-time limit. 



A. Thomas- Fermi parabolic density profiles 

In this subsection we consider specific density profiles 
of the molecular BEC and use the perturbative results of 
Eqs. ([26]) and ([27]) for calculating atom-atom pair cor- 
relations and the relative number squeezing in the short- 
time limit. As one of the most typical situations, we 
start with a parabolic density profile characteristic of an 
interacting BEC in a harmonic trap in the TF limit. 



1. One dimension (ID) 

In ID we assume that the molecular BEC profile is 
given by an inverted TF parabola, po(x) = PoO- — 
x 2 /R' T - F ) for x < Rtf [and po(x) — for x > Rtf], 
where pa = po(0) is the peak density. Accordingly, the 
effective coupling constant is given by 17(3;) = x v / po(l — 
x 2 / R^p) 1 / 2 . The integrals appearing in Eqs. ([2"rJ]) and 
([27]) can be expressed in terms of Bessel functions, using 
the following integral representation [72| : 



2(9/2)" 



+ 1/2) 



/ df(l-£ 
Jo 



2\v-l/2 



cos(gC), (28) 



where v > —1/2 and Y(y) is the gamma function. Ac- 
cordingly, the integrals in Eqs. ([2"u]) and (12~T]) yield 

\mi 2 (k,k ,t)| ~ — : v — , (29) 



(k + k')R 



TF 



rij(k, k' , t) 



V2t 2 x 2 PoR T F J3/2 ((k - k')R TF ) 



[(k - k')R 



TF 



3/2 



(30) 



150 (a! 



(b) CL correlation 

bosons 




Figure 10: (Color online) Analytic results of Eqs. (|3ip and 
(|32[) for the atom-atom pair correlation functions at t — O.lto, 
for the case of a TF parabolic profile of the molecular BEC 
with Rtf = 250 fj.ni. Other parameters are as in Fig. [2] (a) 

<?12 {k, ko, O.lto) as a function of k for fermionic and bosonic 
atom pairs in two different spin states, showing the back-to- 
back correlation peak at k — — ko; (b) pair correlation for the 
same-spin atoms (k, ko, t) in dissociation into distinguish- 
able (fermionic or bosonic) atom pairs, showing a CL HBT 
peak for bosons or a HBT dip for fermions at k — ko- 



Substitution of these expressions into Eqs. ([2"0"]) and 
([21]) leads to the following results for the BB and CL 
pair correlation functions: 



g^(k,k',t) 



9 { 3{k,k',t) 



9tt 2 



[Ji ((* + k')R TF )Y 



16t 2 x 2 p [(k + k')R TF y 



1± 



9tt [h/2{{k~k')R T F)\ 
2 [(k - k')R TF ] 3 



(31) 



(32) 



These results are valid for t <C to, and are plotted in Fig. 
[TU1 as a function of k, for k' = k and t = O.lt . Once 
scaled with respect to the corresponding peak values and 
plotted as in Fig. \7\ the curves in Fig. [TU] follow closely 
the bosonic and fermionic numerical results shown in Fig. 
[7] Interestingly, even beyond the strict range of applica- 
bility of the analytic results of Eqs. ([3"T]) and ([3"2"|) . they 
show good agreement with the numerical results of the 
previous section, valid for up to t ~ to. For example, at 
t = to the differences between the BB and CL correla- 
tion widths, evaluated using the analytic and numerical 
results, are less than 10%. 

The correlation widths that follow from Eqs. ([2"0]) and 
(JUJ) are 



; (bb) 

„(CL) 



w ~ 1.62/J&TF, (33) 
1.1210^1.81/^, (34) 



and therefore W ( CL )/ W ( BB ) ~ 1.12. Here w ~ 1.62/i? T F 
is the width of the molecular BEC momentum distribu- 
tion, n (k) — I J dxyj po{x) exp(— ifcx)/(27r) 1 / 2 | 2 , given 
by 



n (k) 



npo [Ji(kR TF )Y 
2 k 2 



(35) 
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Inverting the relationship between the second- 
and first-order correlation functions, \g^\k,k' ,t)\ — 

\J \9jf{k, k',t) — 1|, we can also find the width of the 

first-order correlation function, Afc( coh ', which gives the 
phase coherence length and defines the size of an atomic 
mode in free space. From the above analytic result 

(2) 

for g\J (fc, fc' ,t), we find that the first-order coherence 
length in the short-time limit is approximately Afc( coh ) = 
2.50/i?TF- In terms of the momentum width of the source 
condensate w, this corresponds to Afc( coh ) = 1.54u>, and 
shows that the first-order coherence length is larger than 
the second-order correlation length w^ CL \ 

Using the asymptotic behavior of the Bessel function, 
J v (z) ~ («/2)"/r(i/+l) for z < 1 (v ± -1, -2, . . .), Eqs! 
(|29| and (|30| at k' = =Ffc give, respectively, the peak value 
of the anomalous density m\2(k 1 t) = m\2(k, — k, t) and 
the atomic momentum distribution in the spin state j, 
rij(k, t) = rij(k, k, t): 



\mi 2 {k,t)\ 
rij(k, t) 



tX\/PoRTF 
4 

2t 2 X 2 paRTF 

3tt ' 



(36) 
(37) 



The peak densities are uniform in the short-time limit, 
corresponding to spontaneous initiation of dissociation, 
which populates the atomic modes uniformly without the 
need to strictly conserve energy. A double-peaked struc- 
ture with maxima around \k\ = ±fco (which in 2D and 3D 
turns into a circle and a sphere of radius |k| = ko [3(|[39j|) 
forms later in time [l7j |. when the present analytic results 
are not applicable. 

Similarly, we can find that the peak values of the pair 
correlation functions (|31|) and (|3"2"T) are given by 



gg\k,-k,t)~ 1 + 



9tt 2 



Q4t 2 x 2 Po ' 



9§\k,k,t) 



2, bosons, 
0, fermions. 



(38) 



(39) 



The inverse square dependence of the BB correlation 

(2) 

peak, g\ 2 '(k,—k,t), on time (for t <C to) is in good 
agreement with the numerical results of the undepleted 
molecular-field approximation, as well as with the re- 
sults of first-principles simulations using the positive- P 
method (see Fig. |6|). The peak values of the CL cor- 

(2) 

relation, g\-'(k,k,t), correspond to the expected HBT 
bunching for bosons and the HBT dip for fermions. 

Finally we calculate the short-time asymptote for the 
relative number variance [Eq. (j24|) ] and find the following 
simple result: 



,(*) 



1 



37rAfci?xF 
32 ' 



(40) 



In this limit the atomic mode populations are much 
smaller than 1, bosonic and fermionic shot noises are ap- 
proximately equal to each other, and we do not see any 



difference in the relative number squeezing for bosons 
and fermions. In Fig. [SJ the results of Eq. (|40p for three 
different sizes of the initial molecular BEC are shown as 
horizontal dashed lines. As we see, the relative number 
squeezing at short times depends merely on the size of the 
source condensate Rtf and the size of the counting cell 
Afc. The small geometric prefactor in the second term is 
determined by the shape of the source condensate (TF 
parabola in the present case) . Together with the resolu- 
tion requirement of Afc < 1/R^f, the smallness of this 
prefactor ensures that Vk ,-k (t) > 0; at the same time 
this implies that the raw (unbinned) squeezing can be 
very weak for small Ptf- 



2, Two dimensions (2D) 

In 2D, the TF density profile of the molecular BEC in 
a harmonic trap is given by po( x ) = Po(l — x 2 /i? TF x — 

V 2 /R\F,y) for X 2 /R%F,x + y 2 /^TF,y < 1 [ an d Po(x) '= 

otherwise], and therefore g(x) = X\/p~o{^ ~ i 2 /^tfi — 
y 2 /i? TF y) 1 ^ 2 ■ The analysis of the short-time asymptotic 
behavior of the correlation functions is essentially the 
same as in ID, except that one has to specify in advance 
the direction of the displacement Ak between the pair of 
momentum vectors k and k' for which the correlations 
are being analyzed. For definiteness, we consider BB 
and CL correlations for which the displacement Ak is 
along one of the Cartesian coordinates, fcj [i = x,y). 
In other words, we consider correlations between k and 
k' = ±k + e, Afci, where is the unit vector in the fcj 
direction; the plus sign is for the CL correlation, and the 
minus sign is for the BB correlation. In the results that 
follow, the dependence on k is absent and we omit it for 
notational simplicity. 

The 2D integrals in Eqs. (f2l))) and (|2"T|) can again be 
performed in terms of Bessel functions (see Appendix |A"| , 
and we obtain 



\mi2(ki,k'i,t)\ 



Uj(ki, fcj, t) 



tXy/PoRTF J3/2 ((ki + ki)RTF,i) 



2V2tt [(fc, + k[)R TF , 



i3/2 



t 2 X 2 p RTF J2 {{h - fcQflxF.i) 

* " [(ki - ki)RrF,i} 2 



(41) 



. (42) 



where Rtf = {Rtf, x Rtf, y) 1 ' is the geometric mean TF 
radius. The preceding expressions lead to the following 
explicit results for the second-order correlation functions: 



g ( $(k i ,k' i ,t)~l + 



Sn [■h/2{(h + k' l )R T¥ , i )\ 
t 2 x 2 Po [(ki + k[)R 



TF,i\ 



1 ± 



64 [J 2 {(kj -kj)RTF,i)Y 

[(ki - fc0i?TF,i] 4 



(43) 



(44) 
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(2) 

The dependence of g\ 2 '(ki, k[, t) on ki + k' { and its peak 
at kj + k[ = implies that it describes the BB correlation 
of different spin atoms with equal but opposite momenta, 
for which k' ~ — k and is offset from — k by an amount 
Aki in the ith direction. Similarly, the dependence of 

(2) 

9jj (^i) k%, t) on ki — k[ and its peak at ki — k[ = cor- 
responds to the CL correlation between pairs of atoms 
in the same spin state, for which k' ~ k and is offset 
by an amount Aki. The qualitative behavior of the BB 
and CL correlation functions is the same as in ID (Fig. 
mjj) whereas the quantitative differences enter through 
the width of the correlations and their respective peak 
values. 

The widths of the BB and CL correlation functions in 
2D are 



w 



■w 



(BB) 
i 

(CL) 



vn ~ 1.81/-R T F,i, (45) 
l.lOwi ~ 1.99/i? TF ,i, (46) 



and therefore w J - CL ' l /?i; J - BB ' 1 ~ 1.10. Relative to the 
source width Wj, the CL correlation width in 2D is nar- 
rower than in ID, while the BB correlation width is 
equal to to,. Here, u>; ~ 1.81/i?TF,i is the width of the 
momentum distribution of the molecular BEC along i, 
found from no(fcj) = \J c? 2 xy / p (x) exp(— ikiXi) / (2ir)\ 2 
(see Appendix as 



n (h) 



i?TF 4 [J3/2(kiRTF,i)] 



(kiR TFti f 



(47) 



The first-order coherence length in 2D, which follows 

33 



from \gf) {kuk'^ = J\gfhk u k[,t) - 1|, is given by 



Ak 



Dh) 



2.75/i? TF ,i and therefore Afc l (coh) ~ 1.52^. 



The peak values of the anomalous and normal densi- 
ties, mi2(k, t) = mi2(k, — k, t) and jij(k, i) = ?ij(k, k, t), 
are given by 



mi 2 (k,t) 



Gtt 



t A x PqRtf 
8tt 



(48) 
(49) 



(2) 

The peak CL correlation function gSS (k, k, t) is the 
same as in Eq. (|39|) . whereas the peak BB correlation is 
given by 



where Afc = (Ak x Aky) 1 / 2 is the geometric mean lattice 

2 

spacing, with Afc giving the counting area. Comparing 
this result with Vfc 0) _/- (£) in ID [Eq. (|4"0"|) ]. we see that 
the raw squeezing is weaker in 2D than in ID, for the 
same size of the molecular BEC and the same lattice 
spacing. 



3. Three dimensions (3D) 



In 3D the TF density profile is given by po( x ) = 
p (l - x 2 /R 2 F ^ x - y 2 /R 2 TF , y - z 2 /R 2 TF J for x 2 /R 2 TF ^ x + 
y 2 /^TF, y + z 2 /^tf, 2 < 1 [ an d Po(x) = otherwise], 
and therefore g(x) = Xy / Po( 1 ~ x2 / r tf,x ~ V 2 l R TF,y ~ 
z 2 jR\ F J 1 ' 2 - As in 2D, we are interested in BB and CL 
density correlations between two momentum components 
at k and k', for which the displacement Ak = k — k' 
is along one of the Cartesian coordinates, ki, where 
i = x,y,z. The 3D integrals in Eqs. ((2"6"|) and (f2"7) can 
again be performed in terms of Bessel functions (see Ap- 
pendix [Bj , and we obtain 



|mi 2 (fci,fci,i)| 



tX^/PoRTF J 2 {{h + fc-)i?TF,i) 



-in 



[(h + kQR 



(52) 



TF,: 



nj (ki, k, t) a * V ^ TF ^ ((fc ; , (53) 



where Rtf — (Rtf.xRtf^Rtf.z) 1 ^ 3 is the geometric 
mean TF radius. The BB and CL correlations follow- 
ing from these expressions are 



(2)/,. ,„/^i, 225n2 [J2{{h + k>)R TF /"'' 



g$(k i ,k! i ,t)~l + 



m 2 x 2 Po [(h + K)r tf A 4 



a m (k y ,w, 1 + 2257r [^((kj-k'^RTF, 



(54) 



(55) 



As in 2D, the qualitative behavior of the BB and CL cor- 
relation functions is the same as in ID (Fig. ITU)) , whereas 
the quantitative differences enter through the width and 
the peak values. 

The widths of the BB and CL correlations in 3D are 
given by 



^(k.-k.t)-! 



1G 



Qt 2 x 2 po 



(50) 



This result is qualitatively similar to the ID result of Eq. 
(|3"8")l . except that the numerical prefactor in the second 
term is different. The same is true for the relative number 
variance. In 2D it is given by 



Vk ,- ko (t)^l 



2(Ak) 2 (R TF ) 2 
9?r 



(51) 



W f B) = Wl ~ 1.99/i?TF, 



w 



j CL) ~ 1.08wi ~ 2.16/^tf,, 



(56) 
(57) 



and therefore w^^'/w^ 3 ' ~ 1.08. The CL corre- 
lation width in 3D relative to the source width tvi 
is smaller than in ID and 2D, whereas the rela- 
tive BB correlation width is the same. Here Wi ~ 
1.99/i?TF,i is the width of the momentum distribution 
of the molecular BEC along i, found from no(kt) = 



15 



| / d 3 x-^/p (x) exp(-ifciXj) / (2n) 3/2 \ 2 (see Appendix |BJ) 



n (ki) 



TTPqRtf [J2(hRTF,i)Y 



(58) 



2 (fc^TF,) 4 

The first-order coherence length in 3D, following 

1|, is given by 



from \gf){hX,t)\ = ^ \ 9 f? {h,k[,t) 



Afc 



(coh) 



2.99/i?TF,i and therefore Afc, 



(coh) 



1.50wi 



This is again larger that the second-order CL correlation 
width wf L) ~ l.OSu;,-. 

The peak values of the anomalous and normal densities 
in 3D are given by 



m 12 (k,t) 



F\ PqRtf 
15tt 2 



^X^/Po^TF 
32tt 



whereas the peak BB correlation is 



fl &'(k,-k,t)~l 



32 2 t 2 x 2 Po' 



The peak CL correlation function gr. . (k, k, t) is the same 

as m Eq. d39]). The result for (k, — k, t) is qualita- 
tively similar to the ID and 2D results of Eqs. |38|) and 
(fSUj) . except that the numerical prefactor in the second 
term is different. 

Finally, the relative number squeezing in 3D is deter- 
mined by 



(59) 
(60) 

(61) 



15(Afc) 3 (fl TF ) 
322 



(62) 



where Afc = (AtAfcyAfcz) 1 ' 3 is the geometric mean lat- 
tice spacing, with Afc 3 giving the counting volume. The 
raw squeezing in 3D is weaker than in ID and 2D, for 
the same size of the molecular BEC and the same lattice 
spacing. 

B. Gaussian density profiles 

In this subsection we present the short-time analytic 
results for the correlation functions in the case of a 
Gaussian density profile of the source molecular BEC, 

po( x ) = Pa exp(— Y^i=i x i/2Sf), where po is the peak 
density and Si is the rms width in the ith direction. 
The results for ID, 2D, and 3D systems can be combined 
through setting D = 1, 2, or 3. As in the case of the TF 
density profile, we are interested in BB and CL density 
correlations between two momentum components at k 
and k', for which the displacement Ak = k k' is along 
one of the Cartesian coordinates, fcj, where i = x, y, z in 
3D (in 2D, i = x,y, while in ID, i = x). We define the 
momentum width of the molecular condensate along fcj 



via <Xj, which corresponds to the rms width 0{ = l/2Si 
of a Gaussian 

n o( k i) = 2 D(^)2D ex P(~ k i / 2cr i). ( 63 ) 

where no (fcj) is defined according to Ho (fcj) = 
\Jd D x^ ~(x)exp(--ihx. l )/(2TT) D / 2 \ 2 and a 



YiiLi a ij i s the geometric mean width. 
With these definitions, the integrals in Eqs. ([2"o| and 
(E3, with g(x) = aVpo(x), give 

txVPo 



l/D 



\m 12 (h,klt)\ ~ 2Z ^_ exp [-(fcj + *0 a /4o?] , (64) 

Ujih, k'^t) ~ 9D + D * 2 ^° D/2 _ exp [-(fcj - fc-) 2 / 8a "i] > 

(65) 

and therefore the BB and CL correlation functions are 
g$(k, k', t) ~ 1 + exp [-(fcj + fcJ) 3 /2o?] , (66) 



4' J (fcj, fci t) ~ 1 ± exp [-(fcj - fc,') 2 /4^ 2 ] • (67) 
The BB and CL correlation widths are 



a( BB ) = a i = 1/25,, 
a {CL) = V2<Ti = l/V2Si 



(68) 
(69) 



The CL correlation width, relative to the source width Oj , 
in the present Gaussian case is broader than in the case 
of a TF parabolic density profile of the source molecular 
BEC, whereas the BB correlation width is equal to the 
width of the source as before. 

The relative number squeezing is given by 

y ko ,_ ko (0 ~ i - (^j D/2 (M) D (S) D , (70) 



where S = ( H2=i &i) is the geometric mean width 
of the molecular BEC in coordinate space and Afc = 



l/D 



YiiLi AfciJ . The general qualitative conclusions 

about relative number squeezing in ID, 2D, and 3D re- 
main the same as for the TF parabolic density profile of 
the source molecular BEC. 



l/D 



VI. EFFECTS OF MOLECULAR DEPLETION 
AND COLLISIONAL INTERACTIONS 

In this section we discuss the role of molecular deple- 
tion and s-wave scattering interactions. We treat these 
effects explicitly using two alternative phase-space rep- 
resentation techniques: first-principles simulations using 
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the positive-P method [73| and a truncated Wigner ap- 
proximation [2I, [75[ . Owing to the fact that these phase- 
space methods are currently well established for bosonic 
fields (see, e.g., [H |H Hf |H, ES M, Iz3 ) > we restrict 
our study only to dissociation into bosonic atoms. De- 
velopment of similar techniques for fermions is underway 
[4ll |78| , but they are so far limited to treating homo- 
geneous systems and therefore are not adopted yet to 
the present problem of dissociation of spatially inhomo- 
geneous molecular condensates. 

To understand the role molecular depletion by itself, 
we first treat the dissociation dynamics governed by the 
Hamiltonian (TTJ) without the s-wave scattering interac- 
tions. Since most of our previous numerical and analytic 
examples were given for the case of distinguishable atoms 
(which has its fermionic counterpart), we examine the 
role of molecular depletion in the same examples. 

In the second part of this section we analyze the role of 
collisional interactions between the atoms and molecules 
and restrict our study to the case of indistinguishable 
bosonic atoms. More specifically, we treat molecule- 
molecule, molecule-atom, and atom-atom s-wave scatter- 
ing interactions described by the respective s-wave scat- 
tering lengths aoo, aoi> and an. Our treatment automat- 
ically incorporates molecular depletion since the molec- 
ular field in all cases is treated quantum mechanically, 
without invoking the mean-field approximation. The re- 
striction to the case of indistinguishable atoms is moti- 
vated by the need to keep the parameter space manage- 
able, while still giving us an overall quantitative under- 
standing of the role of these collisional processes. For 
comparison, in the case of dissociation into distinguish- 
able atoms, a generic treatment would have to incorpo- 
rate six different types of intra- and interspecies scatter- 
ing processes described by six (generally different, and 
yet unknown for most of the species) scattering lengths 
Qij = 0, 1, 2), which is a challenging task and is be- 
yond the scope of the present article. 



A. Role of molecular depletion 

To model the quantum dynamics of dissociation be- 
yond the undepleted molecular-field approximation, we 
use a first-principles phase-space method based on the 
positive-P representation. In this method, each pair of 
the field operators *i(x,i) and *J(x, i) (i = 0,1,2) in 
the Hamiltonian ([T]) in ID is represented by two complex 
stochastic fields *i(x,i) and *i(x,i) whose dynamics is 
governed by a set of stochastic differential equations, 
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(71) 

Here & = (& + i£ 2 )/V2 and £ 2 = (£, + < 4 )/\/2 are 
the complex noise terms, with Q(x,t) (j = 1,2,3,4) be- 
ing real independent Gaussian noises with zero mean, 
(Cj(x,t)) = 0, and the following nonzero correlations: 
(Ci(x,t)Cj(x',t')) = 5 %] 8{x - x')S(t - If). The stochas- 
tic fields *i(x,£) and *i(x, t) are independent of each 
other [*j(x, i) 7^ **(x,t), t > 0], except in the mean, 
(*i(x,i)) = (**(x,i)), where the brackets refer to 
stochastic averages with respect to the positive-P dis- 
tribution function. In numerical realizations, this is rep- 
resented by an ensemble average over a large number 
of stochastic realizations (trajectories) . Observables de- 
scribed by quantum mechanical ensemble averages over 
normally ordered operator products have an exact corre- 
spondence with stochastic averages over the fields *(x, t) 
and *(x,t). 

The initial condition for our simulations is a vacuum 
state for the atomic fields and a coherent state for the 
molecular condensate, with *o(x,0) = *q(x,0). In the 
numerical examples, we assume that the molecular con- 
densates initially have the same density profiles as those 
used in our simulations within the undepleted molecu- 
lar approximation (see Fig. [1]). We note that in the 
undepleted case, these density profiles were assumed to 
originate from the ground-state solution of the Gross- 
Pitaevskii equation in a harmonic trap, which in the TF 
limit gives inverted parabolas whose TF radii depend on 
the strength of molecule-molecule interactions. In the 
present case, we treat the molecular depletion but ignore 
the molecule-molecule s-wave scattering interactions; ac- 
cordingly, the ground state of a harmonic trap would re- 
sult in a Gaussian density profile rather than an inverted 
parabola. To make the present treatment self-consistent, 
yet directly comparable with the numerical examples an- 
alyzed in the undepleted case, we therefore assume that 
the same, near-parabolic initial density profiles are cre- 
ated by tailoring the shape of the trapping potential. 

The results of the exact positive-P simulations are 
shown in Figs. O HI and |H1 the comparison with the re- 
sults of the undepleted molecular approximation is dis- 
cussed in the respective parts of text in Sec. IIV1 A known 
drawback of the positive-P method is that it suffers from 
increasingly large sampling errors due to the boundary 
terms problem [79l [80j as the simulation time increases, 
eventually leading to diverging results. In the present 
examples without the s-wave scattering terms, the useful 
simulation times were limited to t ~ 6.5to. 

In order to go beyond the simulation durations achiev- 
able via the positive-P method, we have also performed 
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simulations using the truncated Wigner-function ap- 
proach. Unlike the positive-P method, it is an approx- 
imate approach as it involves neglecting or truncating 
third- and higher-order derivative terms in the evolution 
equation for the Wigner function. This is necessary in or- 
der to obtain an equation in the form of a Fokker-Planck 
equation which can then be mapped onto a set of stochas- 
tic differential equations. These equations formally ren- 
der as deterministic mean-field (or Gross-Pitaevskii-like) 
equations, 
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however, their stochastic nature and quantum fluctua- 
tions are included by way of a noise contribution in the 
initial state for the molecular and atomic fields. The 
addition of this initial vacuum noise ensures that the ini- 
tial states of *o and *i.2 represent the Wigner function 
of an initial coherent-state BEC for the molecules and 
an initial vacuum state for the atoms, respectively. The 
corresponding stochastic averages with the Wigner dis- 
tribution function correspond to symmetrically ordered 
operator products, so that the calculation of observables 
represented by normally ordered operator products needs 
appropriate symmetrization. 

The results of our simulations using the Wigner func- 
tion method are shown in Figs. [2] and [6] and are in ex- 
cellent agreement with the exact positive-P results, thus 
reinforcing our confidence in the adequacy of this method 
for treating the problem of molecular dissociation. 



B. Role of collisional interactions 

We now turn to the treatment of s-wave scattering 
interactions in the case of dissociation into indistinguish- 
able bosonic atoms. This case is described by the Hamil- 
tonian ©, together with the additional quartic terms, 
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Here Uqq D \ Uq\ D ^ = u\^ D \ and U^ ^ correspond to the 
effective ID coupling constants describing, respectively, 
molecule-molecule, molecule-atom, and atom-atom in- 
teractions that are proportional to the 3D scattering 
lengths aoo, &oi, and an. In the case of a harmonic 
transverse confinement realizing a cigar-shaped ID sys- 
tem, these constants are given by = 2u±au and 
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where uj i is the transverse har- 



The treatment of the s-wave scattering interactions us- 
ing the positive-P method is a challenging task because 
the problem of growing sampling errors becomes more 
severe than before. The useful simulation time with re- 
alistic physical parameters reduces to sub-milliseconds in 
our examples, which is too short to give any new in- 
sights beyond the undepleted molecular approximation 
(see also Ref. [39]). As discussed in Ref. [43, the most 
reliable method in this situation is the truncated Wigner 
approach, which produces stochastic equations that re- 
main stable for much longer simulation durations than 
the positive-P equations. The Wigner equations in the 
present case are 
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monic oscillator frequency [81| 



and the quantum fluctuations are introduced as previ- 
ously via the noise contributions in the initial state for 
the molecular and atomic fields. 

As an example relevant to practice, we treat 87 Rb 2 
molecules as in the experiments of Ref. [l2| and use 
the ID parameters described in [8l|. In Fig. [TT]we plot 
the total number of free dissociated atoms Ni(t) relative 
to their initial number within the molecular condensate 
2^(0). The two solid curves are reference examples, 
corresponding to the undepleted molecular approxima- 
tion (light gray) and the case when the molecular de- 
pletion is included, but all s-wave scattering interactions 
are still absent (blue curve). The dashed and the dash- 
dotted curves are, respectively, for the cases that include 
only molecule-molecule and molecule-atom interactions 
(the case with only atom-atom interactions was studied 
in HU), while the dotted curve corresponds to the inclu- 
sion of all three types of s-wave scattering terms. As we 
see, all these curves agree with the curve for the unde- 
pleted molecular approximation for durations oit/to < 3. 
We therefore conclude that the inclusion of s-wave scat- 
tering interactions reduces the regime of validity of the 
undepleted molecular approximation to durations corre- 
sponding approximately to 5% conversion, which is about 
twice lower than in the absence of these interactions. 

In order to qualitatively understand the behavior of 
the different curves in Fig. [TTJ we note the follow- 
ing two dominant features: (1) the two curves (dashed 
red and dotted magenta) that correspond to having 
nonzero molecule-molecule interactions grow slower than 
the other curves, and (2) all curves with s-wave scattering 
interactions "bend" at around t/to ~ 5 and do not reach 
the maximum conversion efficiency seen in the solid blue 
curve. 

The slower growth of the curves that include the effect 
of molecule-molecule interactions is explained by the fact 
that the molecular condensate in these examples experi- 
ences additional expansion due to the repulsive s-wave 
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Figure 11: (Color online) Fractional total atom number as 
a function of time, N\(t) /2iVo(0), in dissociation into indis- 
tinguishable bosonic atoms. The timescale to = in 
these examples is to — 0.035 s (8l| . The two solid curves 
(light gray and blue) are reference examples, corresponding 
to the undepleted molecular approximation (labeled) and the 
case with molecular depletion but no s-wave scattering in- 
teractions (ciij = nm). The remaining three curves corre- 
spond to simulations in which we include: molecule-molecule 
interactions with ooo = 3 nm (and ooi = an = 0), dashed 
red curve; molecule-atom interactions with aoi = 3 nm (and 
aoo = an = 0), dash-dotted green curve; and all three s-wave 
scattering interactions with aoo = aoi = an = 3 nm, dotted 
magenta curve. In the undepleted case, we assume a near- 
parabolic initial density profile of the molecular BEC, with the 
peak ID density po — 1.8 x 10 7 m _1 and a radius Rtf — 53 
jj,m in the TF limit (originating from a transverse trapping 
potential frequency of oj± /2n = 54 Hz and molecule- molecule 
3D scattering length of aoo = 3 nm). In the cases with de- 
pletion but no s-wave scattering interactions, we assume a 
Gaussian density profile po(x,0) = po exp(— x 2 /2S%) that has 
the same peak ID density and an rms width of S x = 35 p,m 
which is chosen as to closely follow the near-parabolic profile 
for the interacting case in the central part of the cloud. The 
same near-parabolic initial density profile was used in the ex- 
amples with aoo = 3 nm (dashed red and dotted magenta), 
and the same Gaussian profile was used in the example with 
aoi = 3 nm and aoo = an = nm (dash-dotted green). 
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Figure 12: (Color online) (a) Peak molecular density po(0, t) 
as a function of time. The different curves are as in Fig. 1111 
(b) Molecule-molecule local pair correlation in the center of 
the molecular cloud {x = 0,x' = 0,t) as a function of 
time, for the same curves as in (a). 



scattering interactions. Such an expansion is accompa- 
nied by a faster reduction of the molecular peak density 
[see Fig. [T2T a)] compared to the cases with no molecule- 
molecule interactions. Indeed, as we see in Fig. [T2ta) the 
molecular peak density for the dashed red and dotted ma- 
genta curves drops faster than for the other two curves, 
yet the atom number growths for the same cases is the 
slowest in Fig. 1111 Thus, the reduction in the molecular 
density takes place not only because of the conversion to 
atoms, but also because of the expansion of the molec- 
ular condensate. Our estimates show that at t/t Q ~ 3 
the simple expansion alone would reduce the molecular 
density to 0.95 of the original value, which is a larger 
effect than the reduction due to conversion to atoms. In 
terms of an instantaneous effective coupling x(^o(x,t)} 
(interpreted at the level of a time-dependent mean field) , 
which is similar to g(x) = \\/po{x, 0) used previously in 
the undepleted molecular approximation, the additional 
reduction in the molecular density means that the instan- 
taneous atom-molecule coupling and therefore the rate of 
conversion into atoms is reduced more than in the case 
with no expansion of the molecular condensate. Accord- 
ingly, the dashed red and dotted magenta curves in Fig. 
II II have the slowest growth in the atom number. 

The "bending" of the curves at later times can be ex- 
plained by the effect of phase diffusion due to the s-wave 
scattering interactions. The phase diffusion leads to dy- 
namical dephasing of the phase-matching condition for 
efficient conversion, thus suppressing an exponential am- 
plification in the atom number growth. The characteris- 
tic time scale for phase diffusion in the initial stages of 
dissociation can be estimated via 

2tt 

td — .--(id) 777Tm, ' V 75 ) 
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~ u 00 



\Po 



which gives the following results for the three cases shown 
in Fig. QT] tj/to ~ 4.6 for the case with aoi 7^ 0, 
td/to ~ 9.8 for the case with a 00 ^ 0, and i^Ao — 8.7 for 
the case with all s-wave scattering interactions present. 
Here we have ignored the contribution coming from the 
atomic mean field itself (as the atomic density is ini- 
tially negligibly small compared to the molecular peak 
density) and ignored the fact that the mean-field phase 
shifts are spatially dependent and dynamically changing. 
With these remarks in mind, we find that our order-of- 
magnitude estimates of the phase diffusion time are con- 
sistent with the numerical results seen in Fig. [TTJ More- 
over, a model simulation of a uniform system with the 
same s-save scattering interaction terms as in the rel- 
evant three cases of Fig. [TT] reproduce the trend that 
follows from these simple estimates and the order of the 
curves that experience progressively slower phase diffu- 
sion. 

In addition to monitoring the dynamics of the molecu- 
lar peak density, we have analyzed the second-order cor- 
relation function g^(x,x',t) for the molecular field in 
position space in order to understand the deviation of 
the molecular field from the initial coherent state. In Fig. 
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[T27 b) we show the molecule-molecule local pair correla- 

( 2) 

tion in the center of the cloud 9q (0, 0, t) as a function of 
time for the same examples as in Fig. [TST a). We see that 
the strongest deviation from the coherent state value of 

(2) 

<7oo (0)0, t) = 1 occurs for the case with the strongest 
depletion, corresponding to the absence of s-wave scat- 
tering interactions (solid blue curve). In this case, the 
minimum in the peak molecular density — occurring ap- 
proximately at t/to ~ 6.5 — is the smallest and is closer 
to zero, implying that the quantum fluctuations are no 
longer negligible compared to the mean- field part. As 
a signature of this, we see a respective peak in the pair 
correlation function, <7q (0, 0,t) ~ 3, around the same 
time. In contrast to this, in cases with s-wave scattering 
interactions, the molecular depletion is weaker, the peak 
density is still high enough, and therefore the molecular 
field remains closer to the initial coherent state in terms 

( 2) 

of the second-order coherence, with <7q (0, 0, t) ~ 1. 

Finally, in Fig. Q1J] we plot the back-to-back pair cor- 
relation function for the dissociated atoms in momentum 
space and monitor it as a function of time. The dif- 
ferent curves (solid, dashed, dash-dotted, and dotted) 
are for the same numerical examples as in Fig. [TT1 
showing that the s-wave scattering interactions — while 
certainly affecting the quantitative details — have a rela- 
tively smaller effect on atom-atom correlation functions 
at least for dissociation durations up to t/to ~ 4. More- 
over, the numerical results for t/to ^ 1 are hi surpris- 
ingly good agreement with the short-time analytic re- 
sults (shown as squares for the case of the initial TF 
density profile of the molecular BEC and as circles for 
the initial Gaussian profiles), which in the strict sense 
were meant to be applicable only for t/to <C 1. The 
analytic results for the present case of indistinguishable 
bosons follow from Eq. (|22[) . which for k' = —k reads 

as g^{k,-k,t) = 2 + \m(k, -k, t)\ 2 /[n(k, t)] 2 . Using the 
asymptotic expressions in ID given by Eqs. and (f3"7| 
for the case of the TF parabola, we obtain 

9 ^(k,-k, t ),2 + ^-. (76) 

Similarly, using Eqs. (|6"4"|) and (p5)l for the Gaussian den- 
sity profile, we obtain 

g f)(k,-k,t)~2 + ^-. (77) 
X t Po 

Both these expressions have inverse square dependence 
on time t and are shown in Fig. [13] by the curves corre- 
sponding to squares and circles. 

VII. SUMMARY 

We have analyzed the quantum dynamics of dissocia- 
tion of harmonically trapped Bose-Einstein condensates 
of molecular dimers in one, two, and three spatial dimen- 
sions. More specifically, we have initially examined how 
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Figure 13: (Color online) Back-to-back pair correlation for 

dissociated atoms, ffii (&o, —ko,t), as a function of time. The 
different curves (solid blue, dash-dotted green, dashed red, 
and dotted magenta) are as in Fig. 1111 The black squares 
and circles are the results from short-time analytic solutions 
given by Eqs. (|76[) an <l <|77[) . respectively. 



the spatial inhomogeneity of the molecular condensate af- 
fects the conversion dynamics and the atom-atom corre- 
lations in momentum space in the short-time limit. Both 
fermionic and bosonic statistics of the constituent atom 
pairs were considered. Using the undepleted molecular- 
field approximation, we have obtained explicit analytic 
results for the short-time asymptotic behavior of the 
pair correlation functions and the relative atom number 
squeezing. 

For a TF parabolic density profile of the molecular 
BEC, the correlation functions can be expressed in terms 
of Bessel functions and are determined by the momentum 
distribution of the source molecular BEC. The precise 
relationship is different in one, two, and three dimen- 
sions. We have compared the results corresponding to the 
TF density profile with those corresponding to a simple 
Gaussian and found that the correlation widths relative 
to the source width are smaller in the TF case. We have 
also shown that the relative atom number squeezing in 
a given atom counting volume is determined merely by 
the characteristic size of the molecular BEC and by a 
shape-dependent geometric factor and that the squeez- 
ing improves as the size of the BEC is increased. The 
strength of atom-atom correlations for equal but oppo- 
site momenta and the relative number squeezing are the 
strongest in ID systems, where the mode-mixing is a rel- 
atively weaker effect than in higher dimensions. 

The analytic and numerical results in the undepleted 
molecular-field approximation have been compared with 
exact numerical simulations using the positive- P repre- 
sentation in ID, as well as with simulations using the 
truncated Wigner approach, both of which are suitable 
for treating the molecular-field dynamics and its deple- 
tion. The comparison shows that in the absence of any 
s-wave scattering interactions the undepleted molecular 
approximation is valid for dissociation durations corre- 
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sponding to ~ 10% conversion of molecules into free 
atoms. We have also examined the role of the s-wave 
scattering interactions within the truncated Wigner ap- 
proach and found that the phase diffusion due to these 
interactions can significantly modify the conversion dy- 
namics in the long time limit. Nevertheless, for dissoci- 
ation durations corresponding to ~ 5% conversion, one 
can still neglect these interactions and use the undepleted 
molecular approximation. Furthermore, our results indi- 
cate that a possible route to extending the validity of the 
undepleted molecular approximation and minimizing the 
role of collisional interactions is to operate at larger ab- 
solute values of both the atom-molecule coupling x an d 
the dissociation detuning |A|. 

Even though 5% — 10% conversion efficiencies seem 
small, nevertheless they can produce mesoscopic ensem- 
bles of pair-correlated atoms with interesting quantum 
statistics and nontrivial many-body correlations if one 
starts with large-enough molecular condensates, such as 
containing 10 4 — 10 5 molecules in typical ID configura- 
tions. Accordingly, our simple analytic results in the un- 
depleted molecular approximation provide a useful tool 
for obtaining qualitative insights and better than "order- 
of-magnitude" estimates for realistic inhomogeneous sys- 
tems. For a specific system at hand, these results can 
be further refined using the positive-P and truncated 
Wigner methods as demonstrated in the present study. 
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Appendix A: Correlation integrals in 2D 

We first evaluate the integral in the anomalous den- 
sity |mi 2 (k, k', t)\ [Eq. (|26)) ]. in which we take k' = 
— k + e^Afci, so that the integral takes the form 



\m 12 (k x ,k' x ,t)\ = I dxdy 
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where the integration domain A is given by x 2 / R\ F x + 

y 2 / Rtf y < I- For definiteness, we have chosen the di- 
rection i to be the x axis, without the loss of gener- 
ality. Introducing scaled variables, x' = x/Ptf.x and 
y' = y/RTF,y, and transforming to polar coordinates 
x' = r cos 9 and y' = r sin the integral can be rewritten 
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Using the integral representation of the zeroth-order 
Bessel function [73 
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the integral in Eq. (|A2[) can be brought into the following 
form: 



\m 12 (k x ,k' x ,t)\ 
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By making a variable change r = sin0 and using the 
integral representation of the Bessel function of general 
order [72j], 
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with Re p > — 1 and Rc/i > — 1, we finally obtain 

tXy/p6RTF,xRTF,y 



(A5) 



\m 12 (k x ,k' x ,t)\ 



2V2tt 
J3/2 ((k x + k' x )R TF>x ) 



[(kx + k' x )R TF , x 



,3/2 



(A6) 



If the displacement direction i was along y, we would 
obtain the same results except k x — > k y by defining the 
polar coordinates according to x' = r sin 9' and y' = 
rcos9'. Thus, by replacing in \mi 2 (k x ,k x ,t)\ the k x and 
k' x components with ki and k[ (i = x,y) , we arrive at 
the general result of Eq. (dJ). 

To evaluate the integrals in the normal density 
7ij(k,k',i) [Eq. ([27])]. in which we take k' = k + eiAki, 
and therefore 

rij(k x ,k' x ,t) = ^*_f 9 ° / dxdy 
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we follow the same steps as in evaluating \mi 2 (k Xl k' x , t)\. 
In polar coordinates, the preceding integral can be 
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brought into the following form: 
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Introducing r = sin <j> and using Eq. (|A5|) , the integra- 
tion with respect to cf> gives 
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The same result, except k x — > k v , can be obtained if 
the displacement is along y, and we arrive at the general 
result of Eq. (|32"]l. 

Finally, the calculation of the integral in 



n (ki) 



d 2 -x.\/ p {x) exp(-ikiXi) / (2n) 



(A10) 



follows the same pattern as in \mi 2 (ki, k[, t)\ except that 
the relevant term is squared, and one obtains Eq. (|47|) . 



Appendix B: Correlation integrals in 3D 

We first evaluate the 3D integral in the anomalous den- 
sity |mi 2 (k, k', t)\ [Eq. (|26|) ]. in which we take k' = 
— k + e^Afci. For defmiteness, we choose the direction 
i to be along z, so that the integral takes the form 

\m 12 (k z ,k' z ,t)\ = T^TjT / dxdy dz 
A 

« M), fi- W—W--W-) ' (B1) 

\ ^TF.i ^TF.y n TF,z ) 

where the integration domain A is defined by x 2 j R\ F x + 
y 2 1 'R\ F y + Z 2 /R^ F z < 1- Introducing scaled variables, 
x' = x/Rty.x, y' = y/RTF.y, and z' = z/i?TF.z and 
transforming to spherical coordinates x' = r sin # cosy?, 
y' = r sin 9 simp, and z' = rcos#, the integral can be 
rewritten as 



\m 12 (k z , k' z ,t) | = 
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The integral with respect to p is trivial, while the in- 
tegral with respect to 9 is taken using [72j 
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where Rez/ > —1/2. This gives 
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which in turn takes the form of Eq. ()A5|) if we introduce 
r = sin 4>. Accordingly, we obtain 
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For displacements along i = x or i = y directions, one 
can obtain the same result except that k z is replaced by 
fej by an appropriate rotation of the spherical coordinate 
system, so that the final result for \m\ 2 {ki, k[, t)\ takes 
the form of Eq. (|52|) . 

To evaluate the integral in the normal density 
nj(k, k', t) [Eq. ([2"T|]. in which k' = k + e^ki, we follow 
the same steps as for evaluating m 12 (k, k', t)\. In spher- 
ical coordinates, the integral takes the following form: 
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where we have used Eq. (|B3|) . Introducing r = sin <f> and 
using Eq. (|A5|) . we obtain 
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Generalizing this to the arbitrary displacement direction 
i leads to the final result of Eq. (|53|). 
Finally, the calculation of the integral in 



n (ki) 



d 3 K^p (x) exp(-ikiXi)/(2ir) 3/2 



(B8) 
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follows the same pattern as the evaluation of \mi2(ki, k^, t)\, and we obtain Eq. ([551) . 
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